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Abstract. We propose an extension of a standard stochastic individual-based model in popu¬ 
lation dynamics which broadens the range of biological applications. Our primary motivation is 
modelling of immunotherapy of malignant tumours. In this context the different actors, T-cells, 
cytokines or cancer cells, are modelled as single particles (individuals) in the stochastic system. 
The main expansions of the model are distinguishing cancer cells by phenotype and genotype, 
including environment-dependent phenotypic plasticity that does not affect the genotype, tak¬ 
ing into account the effects of therapy and introducing a competition term which lowers the 
reproduction rate of an individual in addition to the usual term that increases its death rate. 
We illustrate the new setup by using it to model various phenomena arising in immunotherapy. 
Our aim is twofold: on the one hand, we show that the interplay of genetic mutations and phe¬ 
notypic switches on different timescales as well as the occurrence of metastability phenomena 
raise new mathematical challenges. On the other hand, we argue why understanding purely 
stochastic events (which cannot be obtained with deterministic models) may help to understand 
the resistance of tumours to therapeutic approaches and may have non-trivial consequences on 
tumour treatment protocols. This is supported through numerical simulations. 
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1. Introduction 

Immunotherapy of cancer received a lot of attention in the medical as well as the mathematical 
modeling communities during the last decades P El El 1313 M- Many different therapeutic 
approaches were developed and tested experimentally. As for the classical therapies such as 
chemo- and radiotherapy, resistance is an important issue also for immunotherapy: although a 
therapy leads to an initial phase of remission, very often a relapse occurs. The main driving 
forces for resistance are considered to be the genotypic and phenotypic heterogeneity of tumors, 
which may be enhanced during therapy, see mm® and references therein. A tumor is a 
complex tissue which evolves in mutual influence with its environment [8]. 

In this article, we consider the example of melanoma (tumors associated to skin cancer) un¬ 
der T-cell therapy. Our work is motivated by the experiments of Landsberg et al. [9], which 
investigate melanoma in mice under adoptive cell transfer (ACT) therapy. This therapeutic 
approach involves the injection of T-cells which recognize a melanocyte-specific antigen and are 
able to kill differentiated types of melanoma cells. The therapy induces an inflammation and 
the melanoma cells react to this environmental change by switching their phenotype, i.e. by 
passing from a differentiated phenotype to a dedifferentiated one (special markers on the cell 
surface disappear). The T-cells recognize the cancerous cells through the markers which are 
downregulated in the dedifferentiated types. Thus, they are not capable of killing these cancer 
cells, and a relapse is often observed. The phenotype switch is enhanced, if pro-inflammatory 
cytokines, called TNF-a (Tumor Necrosis Factor), are present. A second reason for the appear¬ 
ance of a relapse is that the T-cells become exhausted and are not working efficiently anymore. 
This problem was addressed by re-stimulation of the T-cells, but this led only to a delay in 
the occurrence of the relapse. Of course, other immune cells and cytokines are also present. 
However, according to the careful control experiments, their influence can be neglected in the 
context of the phenomena considered here. 

Cell division is not required for switching, and switching is reversible. This means that the 
melanoma cells can recover their initial (differentiated) phenotype [9]. The switch is thus a 
purely phenotypic change which is not induced by a mutation. The state of the tumor is a 
mixture of differentiated and dedifferentiated cells. One possibility to avoid a relapse is to 
inject two types of T-cells (one specific to differentiated cells as above, and the other specific to 
dedifferentiated cells) as suggested also in [9]. 

In this paper, we propose a quantitative mathematical model that can reproduce the phe¬ 
nomena observed in the experiments of [9J, and which allows to simulate different therapy 
protocols, including some where several types of T-cells are used. The model we propose is 
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an extension of the individual-based stochastic models for adaptive dynamics that were intro¬ 
duced in Metz et al. m and developed and analyzed by many authors in recent years (see e.g. 
ini na nans nans in nansi) to the setting of tumor growth under immunotherapy. Such 
models are also referred to as agent-based models or particle systems. Note that we use the 
term individual for T-cells, cytokines or cancer cells, in particular not for patients or mice. 

These stochastic models describe the evolution of interacting cell populations, in which the 
relevant events for each individual (e.g. birth and death) occur randomly. 

It is well known that in the limit of large cell-populations, these models are approximated 
by deterministic kinetic rate models, which are widely used in the modeling of cell populations. 
However, these approximations are inaccurate and fail to account for important phenomena, 
if the numbers of some sub-populations become small. In such situations, random fluctuation 
may become highly significant and completely alter the long-term behavior of the system. For 
example, in a phase of remission during therapy, the cancer and the T-cell populations drop to 
a low level and may die out due to fluctuations. 

A number of (mostly deterministic) models have been proposed that describe the development 
of a tumor under treatment, focusing on different aspects. For example, a deterministic model 
for ACT therapy is presented in [3]. Stochastic approaches were used to understand certain 
aspects of tumor development, for example rate models m or multi-type branching processes; 
see the book by Durrett m or [221 23., [H| • To our knowledge, however, it is a novel feature 
of our models to describe the coevolution of immune- and tumor cells taking into account both 
interactions and phenotypic plasticity. Our models can help understanding the interplay of 
therapy and resistance, in particular in the case of immunotherapy, and may be used to predict 
successful therapy protocols. The main expansions of our models are the following: 


• Two types of transitions are allowed: genotypic mutations and phenotypic switches. 
The characteristic timescale for a mutation to occur can be significantly longer than a 
timescale for epigenetic transitions. 

• Phenotypic changes may be affected by the environment which is not modeled deter¬ 
ministically as in [TBj but as particles undergoing the random dynamics as well. 

• A predator-prey mechanism is included (modeling the interaction of cancer cells and 
immune cells). The defense strategies of the prey are not modeled by different interaction 
rates as in m but by actually modeling the escape strategy, namely switching. 

• A birth-reducing competition term is included which takes account of the fact that 
competition between individuals may also affect their reproduction behavior. 


The paper is structured as follows. In Section [2] we define the model and state the conver¬ 
gence towards a quadratic system of ODEs in the large population limit. In Section [3| we first 
present an example which qualitatively models the therapy carried out in Landsberg et al. |9j 
(Subsection 3.1). We point out a phenomenon of relapse caused by random fluctuations: 


m 


the phase of remission, the typical number of T-cells is small, and random fluctuations can 
cause their extinction, allowing for a growth of the tumor. As a second example (Subsection 


3.2) we study the T-cell therapy with two types of T-cells. In this case an ever richer class of 
possible behavior occurs: either the tumor is cured, or one type of T-cells dies out, or both 
types of T-cells do, or all populations survive for a given time. Subsection T3 presents our 
choice of physiologically reasonable parameters, and the reproduction of experimental observa¬ 
tions together with predictions. In Section [4] we consider the case of rare mutations. Without 
therapy (Subsection |4.1[ ), we first show how to extract an effective Markov chain on the space 
of genotypes, whose transition rates are given by the asymptotic behavior of a faster Markov 
process on the space of phenotypes. In particular, we define a notion of invasion fitness in this 
setting, and describe how we can obtain a generalization of the Polymorphic Evolution Sequence 
introduced in m ■ Second (Subsection |4.2[ ), we study the interplay of mutation and therapy. 
We consider birth-reducing competition between tumor cells and show that the appearance of 
a mutant genotype may be enhanced under treatment. 
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2. The model 

We introduce a general model, which contains three types of actors: 

• Cancer cells: each cell is characterized by a genotype and a phenotype. These cells 
can divide (with or without mutation), die (due to age, competition or therapy) and 
switch their phenotype. We assume that the switch is inherited by the descendants of 
the switched cells. 

• T-cells: each cell can divide, die and produce cytokines. 

• Cytokines: each messenger can vanish and influence the switching of cancer cells. 

The trait space, X, is a finite set of the form 

X = G xVCZUW = {g u .. -,g\g\} x {pi,... ,p\ P \} 

U {zi,.. .,z\ z \j U {u>i, • • • ,W|w|} (2- 1 ) 

where (g,p) G Q x V denotes a cancer cell with genotype g and phenotype p, z G Z a T-cell of 
type z and w G W a cytokine of type w. We write | • | for the number of elements of a set and 
LJ for disjoint unions of sets. 

2.1. General notations and parameters. A population at time t G M+ is represented by the 
measure 

Vt = J2 v t( x ) s *> ( 2 . 2 ) 

x£X 

where x) is the number of individuals of type x at time t and 6 X denotes the Dirac measure at 
x. Here, I\ is a parameter that allows to scale the population size and is usually called carrying- 
capacity of the environment. The dynamics of the population is described by a continuous time 
Markov process, (u/ 4 )t>o, with the following rates: 

Each cancer cell of type ( g,p ) is characterized by 

• natural birth and death rates: b(p) G R+ and d(p) G M+. 

• competition kernels: c(p,p)K ~ 1 G R+ and Cb(p,p)K~ 1 GM+, where the first term in¬ 
creases the death rate and the second term, called birth-reducing competition, lowers 
the birth rate of a cancer cell of phenotype p in presence of a cancer cell of phenotype p. 
If the total birth rate is already at a level 0, then Cb(p,p)K~ 1 G M+ acts as an additional 
death rate. 

• therapy kernel: t(z,p)K~ 1 G R+ additional death rate of a cancer cell of phenotype p 
due to the presence of a T-cell of type z. In addition, £\ d] (z,p) G No cytokines of type 
w are deterministically produced at each killing event. 

• switch kernels: s 9 (p,p ) G M+ and Sw(p,p)K ~ 1 G R+ denote the natural and cytokine- 
induced switch kernels from a cancer cell of type (g,p) to one of type ( g,p ). 

• mutation probability and law: p g G [0,1] denotes the probability that a birth event of a 
cancer cell of genotype g is a mutation. m((g,p), {g,p)) G [0,1] encodes the probability 
that, whenever a mutation occurs, a cancer cell of type ( g,p ) gives birth to a cancer cell 
of type (g,p). By definition m((g,p), (g,p)) = 0 and 'Eg,p m ((9,P), (. 9,P )) = 1- 

Each T-cell of type z is characterized by 

• natural birth and death rates: b(z) G M+ and d(z) G M+. 

• reproduction kernel: b(z,p)K~ 1 G M+ denotes the rate of reproduction of a T-cell with 
trait z in presence of a cancer cell of phenotype p. In addition, £^ od (z,p) G N cytokines 
of type w are deterministically produced at each reproduction event. 

Each cytokine of type w is characterized by 

• natural death rate: d(w) G M+. 

• the molecules are produced when a cancer cell dies due to therapy or a T-cell reproduces. 

Note that the relation between Q and V is encoded in the switch kernels. They specify which 

phenotypes are expressed by a given genotype and influence the proportions of the different 
phenotypes in a (dynamic) environment. 
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Figure [T] provides a graphical representation of the transitions for a population with trait 
space X = {x = (g,p),y = (g,p')} hJ { z x } U {re}, which constitute our model for the ACT 


therapy described in Landsberg et al. [9], see Subsection 3.1 



FIGURE 1. Dynamics of the process (without mutations) modelling the ex¬ 
periments described in [9]. Here, x denotes differentiated melanoma cells, y 
dedifferentiated melanoma cells, z x T-cells and w TNF-a. At each arrow the 
rate for occurrence of the corresponding event is indicated (e.g. birth is illus¬ 
trated with two arrowheads and death with an arrow directed to f). 


2.2. The dynamics. Let A4 K (X) = $xi '■ n £ N,xi,..., x n G A}, i.e. the set of finite 

point measures on X rescaled by K . Then, for each K E N, the dynamics of the A4 K (A)-valued 
Markov process, (v^) t>0l describing the evolution of the population at each time t, can be 
summarized as follows: 

At time t = 0 the population is characterized by a given measure € M K (X). Each 
individual present at time t has several exponential clocks with intensities depending on its 
trait x 6 X and on the current state of the system. We use the shorthand notation v^(p) = 
Yhg&g v i?{g,p) and [_-J± to denote the positive/negative part of the argument. 


(1) Each present cancer cell of trait (g,p) £ Q xP has 

- a clonal reproduction clock with rate (1 — p g ) b(p) — c biPiP) v t' (P) 

Whenever this one rings, an additional cancer cell of the same trait ( g,p ) appears. 

- a mutant reproduction clock with rate p g b{p) — c b{PiP) u t C ip) 

Whenever this one rings, a cancer cell of trait (g,p) appears according to the kernel 
m((g,p), (, g,p )). 

- a natural mortality clock with rate d(p)+J2p^v c (fb P) iy t < (p)+ b(p) — J2pev c b{p > p) u i < (P) 


Whenever this one rings, the cancer cell disappears. 

- a therapy mortality clock of rate J2zez t{ z iP) l 't [ (z). 

Whenever this one rings, the cancer cell disappears and £^ n (z,p) cytokines of type w 
appear according to the weights t(z , p)^ (z). 

- a natural and cytokine-induced switch clock with rate ^2p^-p{s 9 {p,p)+^2 w ^y^ s^,{p,p)i/^(w)). 
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Whenever this one rings, this cancer cell disappears and a new cancer cell of trait (g,p) 
appears according to the weights s 9 (p,p ) + P2weW s w(p,p)^P ( w). 

(2) Each present T-cell of trait z G Z has 

- a natural birth clock with rate b(z ) 

- a natural mortality clock with rate d(z) 

- a reproduction clock with rate E p ep b(z,p)vp (p). 

Whenever the reproduction clock of a T-cell rings, a additional T-cell with the same 
trait z and lw° d (z,p) cytokines of type w appear according to the weights b(z,p)up (p). 


(3) Each present cytokine has 

- a mortality clock with rate d(w). 

Moreover, an amount of £^ n (z,p) particles (of trait w) are produced every time a T-cell 
of trait z kills a cancer cell of phenotype p, and a number of £w° d (z,p) particles are 
produced every time a T-cell of trait z is produced in the presence of a cancer cell of 
phenotype p. 


The measure-valued process (v/')t>o is a Markov process whose law is characterized by its 
infinitesimal generator L K which captures the dynamics described above (cf. {25 | Chapter 11 
and EDI). The generator acts on bounded measurable functions from Ai K into M, for all 
rj G M k by 


(L K 4>) (r?) 


22 {p + 

(g,p)&gxV 


“(fl.p) 

K 


(! - Pg) 


Kp) ~^2c b (p,p)r](p) 


•p&V 


Krjig.p) 

+ 


+ 22 ipt>(p- % i ) - < ^( r ?)) 


+ E E 

(g,p)£GxV z&Z 


x ( d(p) + 22 KP-/P) r l(p) + 
V p&V 

6 , 


Kp) - 22 c - b (pi p) r i(p) 

p&V 


I Kv(g,p) 


V - "-pf- + E™eW tw l \ z ,P)l? ) - <t>(v) ) t(z,p)ri(z)Krj(g,p) 


+ E E i* ( 1 +‘-V - V) - m) 

(g,p)eGxV peV 

x ( s 9 (p,p ) + Et K v(g,p) 

+ X] + K + E„ewC od (^P)|) - <K^)) ( b(z,p)v(p)\ Kti(z) 

zez pev V / 

+ X] fa + a) “ 0fa)) KP K vK) + £>fa-£)- ^( 77 )) d(z)Krj(z) 

z£Z z^LZ, 

+ 22 Kl ~ 6 it) ~ <X r ?)) d(w)Krj(w) 

w£\V 

+ 22 22 ( ?? +%^) _ ^fa)) 

(gv)&GxV (g,p)eQxV 


x p g m((g,p), (g,p)) 


Kp) ~ 22 Cb(p,p)v{p) 

p'£V 


Kg(g,p)- 

+ 


(2.3) 
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All the figures of this article containing simulations of the above model were generated by a 
computer program written by Boris Prochnau. The pseudo-code can be found in the Appen¬ 
dix El 

2.3. Large population approximation. There is a natural way to jointly construct the pro¬ 
cesses (p t A )t>o for different values of K, see [25] . For this joint realization, the sequence of 
rescaled processes ((v^)t>o)K converges almost surely as K tends to infinity to the solution of a 
quadratic system of ODEs, as stated in the following proposition, which can be seen as a law of 
large numbers. The deterministic system, which provides (partial) information on the stochastic 
system, consists of a logistic part, a predator-prey relation between T-cells and cancer cells, a 
mutation and a switch part. 

Note that the population can be represented as a vector Vx(t) := {y^{x)) x& x of dimension 
d = \G\-\P\ + \Z\ + \W\. 

Proposition 1. Suppose that the initial conditions converge almost surely to a deterministic 
limit, i.e. limx^-oo Vk(0) = u(0). Then, for each T E K+ the sequence of rescaled processes 
(yi<{t))o<t<T converges almost surely as K —> oo to the d-dimensional deterministic process 
which is the unique solution to the following quadratic dynamical system: 

For all ( g,p ) E Q x V , 


dn 


(g.p) 

dt 


n (s,p) 1 (1 hg) 


b (P ) - Y c b(P’P) n (g,p) 
{g,p)£Gx'P 


KP)~ Y c b(P’P) n (g,p) 

{gp)&QxV 


d (p) ~ Y c (P’P) n (~9,p) 

{g,p)eGxV 


Y K^p) n z -Y\ s9 (P’P) + Y s w(P’P) n u 


z£Z 


P&V 


ioEW 


+ Y n (9,p) s3 (p,p)+ Y 


p&V 


wGW 


+ Y n r9,p)[p9 m ii^p)^9,p)) 

(9iP)sSx"P V 


b (p)~ Y c b(PiP') n {g',p') 

{g',p')^QxV 


for all z E Z, 


for all w E W, 
dn„, 


dtL 

dt 


n Z lb(z)-d(z) + Y b ( z ’P) n (g, P ) )> 

V [g,p)&QxV / 


(2.4) 


dt 


= - n, 


4{w) + y n (g,p) Y (fi ll i z ip) l i z 'P) +c° d (^»p) K z >p)) n z- 

(g,p)£Gx'P 


More precisely, P (lim^-^oo sup 0<t<T | Vx(t) — n(t) | = 0) =1. where n(t) denotes the solution to 
Equations \2A\ with initial condition u(0). 

Proof. This result follows from Theorem 2.1 in Chapter 11 of [25]. □ 


It is an important feature of stochastic models opposed to deterministic ones that populations 
can die out. There are two main reasons for the extinction of a population for finite K: first, the 
trajectory of the population size is transient and passes typically through a low minimum. In this 
case, random fluctuations can lead to extinction before the population reaches its equilibrium. 
Second, fluctuations around a finite equilibrium cause extinction of a population after a long 
enough time. The second case happens at much longer times scales than the first one. In both 
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cases, the value of K plays a crucial role, since it determines the amplitude of the fluctuations 
and thus the probability of extinction. 

The relevant mutations in the setup of cancer therapy are driver mutations and appear only 
rarely. In this case, more precisely, when the mutation probabilities, p g = /j^", tend to zero as 
I\ —> oo, mutations are invisible in the deterministic limit. Due to the presence of the switches 
the analysis of the system is difficult. It is not a generalized Lotka-Volterra system of the form 
n = n/(n), where / is linear in n. In Section [4] we show how to deal mathematically with rare 
mutations and their interactions with fast phenotypic switches or therapy. 


3. Relapse due to random fluctuations 


3.1. Therapy with T-cells with one specificity. In this section we present an example 
which qualitatively models the experiment of Landsberg et al. [9], where melanoma escape ACT 
therapy by phenotypic plasticity in presence of TNF-a. Mutations are not considered here, 
since this was not investigated in the experiments. 

Let us denote by x := ( g,p ) the differentiated cancer cells, by y := (g,p') the dedifferentiated 
cancer cells, by z x the T-cells attacking only cells of type x, and by w the TNF-a proteins. 
We start with describing the deterministic system and denote by n,; its solution for trait i G 
{x, y, z x , re}. It is the solution to the following system of four differential equations: 

n x = n X (b(x) - d{x) - c(x, x)n x - c(x, y)n y -s{x,y ) - s w (x, y)n w - t(z x , x)n Zx ) + s(y,x) n y 


(3.1) 


n y = 

n y 

{Ky) - d{v) - c (y 

, y)n y - C 

(y,x)n x - 

s(y,x)) + s(x,y) n x 

+ s w (x, 

y )' 

*z x = 

n 2; 

c (b(z x ,x)n x - d(z 3 

0) 





— 

-n W d(w) + (£^(z x , 

x)t(z x ,x 

) + C od (^ 

xi x) b(z x , x')')n x n Zx . 



We fix the 

following parameters: 







b(x) = 3 

Kv) = 

3 

b(z x , x) = 8 

/'kill ( 7 

x) 



d{x) = 1 

d{y) = 

1 

t{z x ,x ) = 28 

t'W \^x 




c(x, x) = 0.3 

c{y,x ) 

= 0 

d(z x ) = 3 

d(w) = 

15 



c(x,y) = 0 

c(y,y ) 

= 0.3 


s w (x,y) 

= 



s{x,y) = 0.1 

s{y,x) 

= 1 




and initial 

conditions: 









(n x ,rij / ,n^,n lu )(0) 

= (2,0,0,0.05,0). 




= 0 


(3.2) 


(3.3) 


The solution to the deterministic system <ED with parameters (|3.2[) and initial conditions (|3.3[) 


can be seen on Figure [4] (A). There are three fixed points in this example, see the black dots 
on Figure [3] (B): Poooo where all populations sizes are zero, P xy 00 where the T-cells and TNF-a 
are absent and both melanoma populations are present and P xyZx w where all populations are 
present. P xyZx w is the only stable fixed point and P xy 00 is stable in the invariant subspace 
{n Zx = 0} (i.e. if the T-cell population is zero). To highlight the qualitative features of the 
system, we choose parameters such that the minimum of the T-cell population during remission 
is low, and such that the equilibrium value of melanoma of type x in presence of T-cells is low, 
whereas equilibrium values of both melanoma types in absence of T-cells are high. 

For initial conditions such that the number of differentiated melanoma cells, n x (0), is large, 
the number of injected T-cells, n 2a ,(0), is small, and the numbers of dedifferentiated melanoma 
cells, %(0), and TNF-a molecules, 1x^(0), are small or equal to zero, the deterministic system 


is attracted to P 


x y z x w 


the T-cell population, n Zx , increases in presence of its target x , TNF-a 
is secreted, and the population of differentiated melanoma cells, n x , shrinks due to killing and 
TNF-a induced switching, whereas the population of dedifferentiated melanoma cells, n y , grows. 

For the stochastic system, several types of behavior can occur with certain probabilities: 
either the trajectory stays close to that of the deterministic system and the system reaches a 
neighborhood of the fixed point P xyZxW (Fig. [ 2 ] (A)) or the T-cell population, u K (z x ). dies out 
and the system reaches a neighborhood of P xy 00 (Figure [ 2 ] (B)). In the latter case the TNF-a 
population, v K (w), becomes extinct shortly after the extinction of the T-cells, u K (z x ), and 
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the population of differentiated melanoma cells, v K (x), can grow again. Moreover, TNF-a 
inducing the switch from x to y vanishes and we observe a relapse which consists mainly of 
differentiated cells. Depending on the choice of parameters (in particular switching, therapy or 
cross-competition), a variety of different behavior is possible. 




FIGURE 2. Simulations of the evolution of melanoma under T-cell therapy for 
parameters (3.21 and initial conditions (3.31. The graphs show the number of 
individuals divided by 200 versus time. Two scenarios are possible for therapy 
with T-cells of one specificity: (A) T-cells z x survive and the system is attracted 
to Px y z x w , or (B) T-cells z x die out and the system is attracted to P xy oo- 


3.2. Therapy with T-cells of two specificities. A therapy can only be called successful if 
the whole tumor is eradicated or kept small for a long time. A natural idea is thus to inject 
two types of T-cells in future therapies as suggested in [Dj. To model this scenario, we just need 



FIGURE 3. Qualitative description of the role of the fixed points. (A) Sketch 
of the invariant subspaces, stability of the fixed points, and schematic repre¬ 
sentation of the dynamics of the deterministic and the stochastic processes. 
(B) Projection of the fixed points onto n x and n y . The green area contains the 
fixed points which correspond to cure or remission and the red area contains 
those describing relapses. 
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to add T-cells attacking the dedifferentiated cells as new actors to the setting described above. 
We denote them by z y . The system contains one more predator-prey term between y and z y : 

n x = n X (b(x) - d{x) - c(x,x)n x -c(x,y)n y -s w (x,y)n w -s(x,y ) - t(z x , x)n z ^j + s(y,x)n y 

n y = ny (b(y)- d(y)~ c(y,y)n y - c(y,x)n x - s(y,x) - t{z y ,y)n Zy ^j + (s w (x,y)n w + s(x,y))n x 
^z x — R (b(z x , x'jxij: d{z x )) 

n Zy = n Zy(b(z y , y)n y — d(z y )) 

n w = -n W d(w) + (l^ n (z x , x) t(z x , x) + £V rod (z x , x) b(z x , x))n x n Zx 

+ (^fta/> v) t{z y , y) + C od (%> y) K*y> y))^z y (3.4) 

(3.5) 

(3.6) 


In addition to parameters (3.2), we use the following ones: 

t(zy,y) = 28 e^ l \z y ,y) = 1 d(z y ) = 3 


b(z y , y) = 14 C° d (%, y) = 0 


and initial conditions: 


(n x ,n y ,n Zx ,n z ,n w )(0) = (2,0,0,0.05,0.2,0). 


The solution to the deterministic system ( |3.4[ ) with parameters (3.2) and (3.5) and initial 
conditions (3.6) can be seen on Figure [4] (B). Initial conditions with only T-cells of type z y or 
without any T-cells can be seen on Figure [4] (C) and (D). 





FIGURE 4. Solutions to the deterministic system (3.4) with parameters (3.2) 
and (3.5). With n(0) := (n x , n y , n Zx , n z ,n u ,)(0) as the initial conditions, the 


system is attracted to the indicated fixed point: 

(A) n(0) = (2,0,0.05,0,0), fixed point P xyZx 0w, 

(B) n(0) = (2,0,0.05,0.2,0), fixed point P xyZxZyW , 

(C) n(0) = (2,0,0, 0.2,0), fixed point P xy0Zy w, 

(D) n(0) = (2,0,0,0,0), fixed point P xy ooo- 
































A STOCHASTIC INDIVIDUAL-BASED MODEL FOR IMMUNOTHERAPY OF CANCER 


11 





FIGURE 5. Simulations of the stochastic evolution of melanoma under T- 
cell therapy for parameters (3.3) and (3.51. The graphs show the number of 
individuals divided by 200 versus time. Possible scenarios for therapy with 
T-cells of two specificities: (A) T-cells z x and z y survive and the system stays 
close to P X yz x z y w, (B) T-cells z x die out and the system is attracted to P xy oz v w, 
(C) T-cells z y die out and the system is attracted to P xyZj; ow, (D) Both T-cell 
types z x and z y die out and the system is attracted to P xy ooo- (E) Transition 
between cases (A) and (C). (F) the tumor is eradicated (corresponding to 
Fooooo)- 


The introduction of z y adds two new fixed points (see the blue dots on Figure[3](B)): P xyZxZyW 
is the new stable fixed point with all non-zero populations, and P xy oz y w corresponds to the 
absence of the T-cell population of type z x . The invariant subspaces are now {n Zx = 0}, in 
which P xy oz y w is stable, {n Zy = 0}, in which P xyZx0w is stable and {n Zx = 0} D {n Zy = 0}, in 
which P xy ooo is stable. Note that P xyZx Ow, corresponding to P xyZx w from the last section, is 
unstable in the enlarged space. 
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With the same initial conditions as before and n^O) small but positive, the deterministic 
system is attracted to the stable fixed point P X yz x z y w : the T-cell population, n Zx , increases in 
presence of its target x, TNF-a is secreted, and the differentiated melanoma population shrinks 
due to killing and switching, the population of dedifferentiated melanoma grows, but is regulated 
and kept at a low level by the T-cells of type z y . Similarly, is regulated by n Zx . 

We choose the parameters such that the minima of the two types of T-cells during remission 
are low, so that they have a large enough probability to die out in the stochastic system. Since 
at the beginning of therapy no or only very few dedifferentiated melanoma cells are present, 
the population of T-cells of type z y starts growing only later. In order to avoid their early 
extinction a higher initial amount of these T-cells can be injected. There are now five main 
different scenarios in the stochastic system (see Figure [5]). Either the T-cells of type z x (B), 
or the T-cells of type z y (C), or both of them die out (D). Also all populations can survive for 
some time fluctuating around their joint equilibrium (A). The fifth scenario is a cure, i.e. the 
extinction of the entire tumor due to the simultaneous attack of the two different T-cell types 
(F). T-cells and TNF-a vanish since they are not produced any more in the absence of their 
target. Of course, transitions between the different scenarios are also possible, e.g. the system 
could pass from Case (A) to (B) or (C) and then to (D), see Figure [5] (E). Furthermore, note 
that setting the switch from x to y to zero introduces an additional scenario: it is then possible 
that a relapse appears, which consists only of differentiated melanoma cells. 

Starting from our choice of initial conditions, the deterministic system converges to P xy z x z y , 
but the stochastic system can hit one of the invariant hyperplanes due to fluctuations, and is 
driven to different possible fixed points, see Figure [3] (A). The transitions between the different 
scenarios can be seen as a metastability phenomenon. 


3.3. Reproduction of experimental observations and predictions. 


3.3.1. Comparison of experimental observations and simulations. The parameters of Section [3] 
are chosen ad hoc to highlight the influence of randomness and the possible behavior of the 
system. Let us now show that our models are capable to reproduce the experimental data 
of Landsberg et al. [9] quantitatively. The choice of parameters is explained below (Subsec¬ 
tion 3.3. 2\. 

Figure |6| (A) shows the experimental data of [9] whereas Figure[6](B) shows the results of our 
simulations. Each curve describes the evolution of the diameter of the tumor over time. In the 
stochastic system two situations can occur: first, the relapse consists mainly of differentiated 
melanoma cells and the tumor reaches its original size again after 90 days. This is the case if 
the T-cells die out. Second, the relapse consists mainly of dedifferentiated cells and the tumor 
reaches its original size again after roughly 190 days. This is the case if the T-cells survive the 
phase of remission, become active again and kill differentiated cancer cells. In the simulations 
the therapy with one type of T-cells pushes the tumor down to a microscopic level for 50 to 
60 days, as in the experimental data. The curves marked ACT in the experimental data in 
Figure [6] (A) are matched by simulation data when the T-cells die out (Differentiated relapse 
in Figure [6] (B)). In the experiments there might be T-cells which lose their function, e.g. due 
to exhaustion, and cannot kill the differentiated melanoma cells. This effect is to be seen as 
included in the death rate of T-cells in the model. They can be re-stimulated and become active 
again which is marked as ACT+Re in Figure [6] (A). Although our model does not include re¬ 
stimulation, the case of surviving T-cells in the simulations (Dediff. relapse in Figure[6] (B)) can 
qualitatively be interpreted as the case of ACT+Re. Note that the scales of the axes are the 
same in both figures and that the experimental findings are met very well by the simulations. 
The simulated curves under treatment start at the beginning of the treatment and not at day 
zero. The detailed pictures showing the evolution of melanoma and T-cell populations during 
the therapy are given in Figure [7j 

As there is no data for the case of two T-cells, numerical simulations of such a therapy 


strategy should be seen as predictions. For the new T-cell population (of type z y ) we choose the 
same parameters as for the first population (of type z x ), just the target is different. The therapy 
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FIGURE 6. Comparison of experimental data obtained by Landsberg et al. 
with simulations for biologically reasonable parameters. The graphs show the 
diameter of the tumor measured in millimeters versus time in days after tumor 
initiation: (A) experimental data, (B) simulated data (K = 10 5 and n Z;r (0) = 
0 . 02 ). 



FIGURE 7. Simulations for biological parameters. Therapy with one T-cell 
type: (A) differentiated relapse (T-cells z x die out), (B) dedifferentiated relapse 
(T-cells z x survive), Therapy with two T-cell types: (C) cure, (D) differentiated 
relapse (both T-cell types die out). 


seems to be very promising: almost all simulations show a cure for these parameters, only very 
few times a relapse occurs. Nevertheless, the behavior of the system (e.g. the probability to end 
up in the different scenarios) depends strongly on the choice of certain parameters, as pointed 
out in the last two sections. In order to give a reliable prediction we need data to obtain safer 
estimates for the most important parameters, which seem to be the switching and therapy rates 
as well as initial values. 

The initial values play an important role for the success of a therapy. In the case of therapy 
with T-cells of one specificity, increasing the initial amount of T-cells has the following effect: 
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the melanoma cells are killed faster, the population of differentiated melanoma cells reaches a 
lower minimum and as a consequence the T-cells pass through a lower and broader minimum. 
The probability that the T-cells die out increases, and a differentiated relapse is more likely 
than in the case of a smaller initial T-cell population. Moreover, the broadening of the minima 
causes a “delay” and both kind of relapses (consisting mainly of differentiated or dedifferentiated 
cells) appear later. But since the extinction of T-cells is more likely, the tumor may reach its 
original size earlier, see Figure [8| For an initial value ten times as large as in Figure [6] (B) the 
probability of an eradication of the tumor is still very small. If the number of T-cells initially 
is half the number of tumor cells, the probability of a favorable outcome is much higher. But 
such a high amount of T-cells is unrealistic. 



FIGURE 8 . Simulations for different initial doses of T-cells: n-^O) = 0.2 and 
n z j0) = 0.02. 


3.3.2. Physiologically reasonable parameters. We explain here how we choose the biological pa¬ 
rameters. Some parameters can be estimated from the experimental data. Recall that the 
subject of [9] is to investigate the behavior of melanoma under T-cell therapy in mice. Without 
therapy the tumor undergoes only natural birth, death and switch events. 

• Choice of birth and death rates: We assume that the number of cells in the tumor is 
described by 

N t ~ No exp(rt), (3.7) 

where iV) denotes the number of cells at time t, Nq the initial population size and r the 
overall growth rate. Note that the estimate of the growth rate is independent of the 
initial value. Figure 4 (A) in the main part shows that the tumor needs roughly 50 days 
(without therapy) to grow from 2 mm diameter to 10 mm diameter. Since the structure 
of a melanoma is 3-dimensional this corresponds roughly to N$o = 1251Vo which implies 
r = 0.1. Unfortunately, no data that allow to estimate the ratio of birth and death 
events are provided. As long as mutations are not considered this should not have a 
big impact and we chose b = 0.12 and d = 0.02 for the differentiated as well as the 
dedifferentiated cells. Landsberg et al. observed that the growth kinetics appear to be 
the same for both cell types, see Supplementary Figure 11 in [9]. 

• Choice of the competition: We assume that the competition has a very little effect here 
because the tumor grows exponentially in the observed time frame and does not come 
close to its equilibrium. We choose the competition between melanoma cells of the same 
type as c(x, x) = c(y, y ) = 0.00005 and between different types of melanoma cells as 
c(x, y) = c(y, x ) = 0.00002. The values are not set to 0 since the melanoma can grow 
only up to a finite size. 
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Choice of the switch parameters: We can now estimate the switching parameters by 
using the data of Supplementary Figure 9e in [9]. In this experiment where cell division 
is inhibited, we can set 6 = 0. Furthermore, the amount of TNF-a is constant and we 
set here n w = 2. Thus, the dynamics of the melanoma populations is described by 

h x =n x (- d(x) - c(x, x)n x - c(x,y)n y - 2 s w (x,y) - s(x,y)) + s(y,x)n y . . 

n y =n y (- d(y) - c(y, y)n y - c(y, x)n x - s(y , x)) + (2 s w (x, y) + s(x, y)) n x 

At the beginning of their observations the switch is very slow and speeds up after 
the first 24 hours. We assume that there is a delay until the reaction really starts and 
thus we choose the proportions at day 1 (n x = 0.81 and n y = 0.19) as initial data and 
choose switching parameters such that roughly the concentrations at day 2 (n T = 0.45 
and n y = 0.54) and 3 (n x = 0.24 and n y = 0.72) are reached as shown in Figure [9j 
Thereby we obtain s(x,y) = 0.0008, s(y, x) = 0.065 and s w (x,y ) = 0.33. Note that 
the experiments we refer to provide only in vitro data and it is not clear if the in vivo 
situation is similar. 



FIGURE 9. Switch in the in vitro experiments for inhibited cell division and 
constant concentration of TNF-a. Dashed lines indicate experimental data. 


Choice of parameters concerning T-cells: It remains to characterize the T-cells. Their 
natural birth rate is set to 0 since they are transferred by adoptive cell transfer and 
not produced by the mice themselves and do not proliferate in absence of targets. We 
assume that they have a relatively high birth rate depending on the amount of cancer 
cells present, b(z x ,x) = b(z y , y) = 2 and produce one TNF-a molecule when they di¬ 
vide, $w oA (z x , x) = $w° d {z y ,y) = 1. Furthermore, we assume that 4.5 cancer cells can 
be killed per hour (including indirect mechanisms), t(z x ,x ) = t(z y ,y) = 108. The rate 
of death for the T-cell population is chosen as d(z x ) = d(z y ) = 0.12. These parameters 
are chosen such that the qualitative behavior of the tumor was recovered. We choose 
the same parameters for the second T-cell type as for the first one because there are no 
data concerning the second T-cell type. 

Choice of starting values and the scale K: We set K = 10 5 , the initial value for the 
differentiated melanoma cell population to 1 and to 0 for the population of dediffer¬ 
entiated melanoma cells. The ratio of differentiated and dedifferentiated cells is not 
known for small tumors which do not result from cell transfer of cells of in vitro cell 
lines. The initial value of the T-cell population is set to 0.02. We assume that the T- 
cells appear directly in the tumor, i.e. the migration phase into the tumor is not modeled. 
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To sum up, biological rates (per day) and initial conditions (in 100 000 cells) are: 


b(x) = 0.12 
d{x) = 0.02 
c(x,x ) = 5 • 10 -5 
c{y,x) = 2 • 10~ 5 
s(x, y) = 0.0008 

Mo) = 1 


b(y) = 0.12 
d{y) = 0.02 
c(x, y) = 2 ■ 10 —5 
cly,y) = 5 ■ 10- 5 
s(y, x) = 0.065 
n a (0) = 0 


b(z x ,x) = 2 
t(z x , x) = 108 
d(z x ) = 0.12 


n % (0) = 0.02 


C T ° d (z x ,x) = 1 

e^ n (z x ,x) = o 

d(w) = 0.2 
s w (x, y) = 0.33 

K = 10 5 


(3.9) 


The additional parameters in the case where a second T-cell is used are: 

t(z y ,y) = 108 £w° d (z y ,y) = 1 d(z y ) = 0.12 


b(z y ,y) = 2 e^ l \z y ,y) = 0 


,( 0 ) = 0.02 


(3.10) 


4. Arrival of a mutant 


The model we defined in Section [2] can be seen as a generalization of one of the standard mod¬ 
els of adaptive dynamics, usually called BPDL, which was introduced in publications of Bolker, 
Pacala, Dieckmann and Law mmm- In a possibly continuous trait space L Cl 11 , the BPDL 
model allows for each individual to reproduce, with or without mutation, or die due to natural 
death or to competition (e.g. this amounts to only keep parameters b(x),d(x),c(x,y),m(x,y) 


and y in Subsection 2.1). 


In the limit of large populations (K — > oo and y fixed) Fournier and Meleard [30j proved 
a law of large numbers: the process converges to a system of deterministic integro-differential 
equations. In the case y = 0 the process converges to the solution of a system of coupled logistic 
equations (of Lotka-Volterra type) without mutations. The limit of large populations ( I\ -A oo) 
followed by the limit of rare mutations (y -A 0) on the timescale t ~ log(l / y) was studied by 
Bovier and Wang m and a deterministic jump process is obtained as limiting behavior. 

The simultaneous limits of large populations and rare mutations, where ( K , y ) -a (oo,0) at 
a rate such that l/(Ky) log/i and a timescale t ~ 1 /(Ky) was studied by Champagnat 
and Meleard mm- At this scale the system has time to equilibrate between two successive 
mutations. The long-term behavior of the population can be described as a Markov jump 
process along a sequence of equilibria of, in general, polymorphic populations. An important 
(and in some sense generic) special case occurs when the mutant population fixates while resident 
population dies out in each step. The corresponding jump process is called the Trait Substitution 
Sequence (TSS) in adaptive dynamics. Champagnat (TH] derived criteria in the context of 
individual-based models under which convergence to the TSS can be proven. The general 
process is called the Polymorphic Evolution Sequence (PES). It is described partly implicitly in 
ED, as it involves the identification of attractive fixed points in a sequence of Lotka-Volterra 
equations that are in general not tractable analytically. 

In situations when the population converges to a TSS, one may take a further limit of small 
mutation steps, a, and look at the time scale t ~ 1/a 2 . The mutant is then of trait y = x + ah 
where h is chosen according to a probability kernel on X . For birth and death rates that vary 
smoothly on A, this ensures a vanishing evolutionary advantage of mutants when a -A 0. The 
TSS converges in this limit to the so-called Canonical Equation of Adaptive Dynamics (CEAD), 
see e.g. [27], which describes the continuous evolution of a monomorphic population in a fitness 
landscape. The convergence of the individual-based model to the CEAD in a single step has 
recently be shown by Baar et al. [253, i- e - limit ( K , y, a) -A (oo, 0,0) of large populations, rare 
mutations and small mutation steps are taken simultaneously, provided 1 /(Ky) log(A')/cr 
and 1/VK C < 1. The time-scale on which this convergence takes place is t ~ l/(I\ya 2 ), 
and corresponds to the combination of the previous two. 

Costa et al. |19| study an extension of the model with a predator-prey relation. The predator- 
prey kernel is an explicit function of parameters describing defense strategies for preys, together 
with the ability of predators to circumvent the defense mechanism. In the limit 1 /( Ky) log K , 
convergence to a Markov jump process generalizing the Polymorphic Evolution Sequence is 
derived, and in the subsequent limit a -A 0, in the case of a monomorphic prey and predator 
populations, convergence to an extended version of the CEAD is obtained. 
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4.1. Rare mutations and fast switches. In this section, we give a generalization of the 
Polymorphic Evolution Sequence in the case of fast switches in the phenotypic space, without 
therapy (no predator-prey term). We consider the case of rare mutations in large populations 
on a timescale such that a population reaches equilibrium before a new mutant appears: 

VV > 0, exp(—I 7K) << /u 2 , as K —>■ oo. (4.1) 

y A log A 


Champagnat and Meleard’s proof of convergence to the PES m relies on a precise study of 
the way a mutant population fixates, which we now describe. A crucial assumption is that the 
large population limit is a competitive Lotka-Volterra system with a unique stable[]] fixed point 
h. The main task is to study the invasion of a mutant that has just appeared in a population 
close to equilibrium. The invasion can be divided into three steps. First, as long as the mutant 
population size is smaller than Ke for a fixed small e > 0, the resident population stays close to 
its equilibrium. Therefore the mutant population can be approximated by a binary branching 
process. Second, once the mutant population reaches the level Ke , the whole system is close to 
the solution of the corresponding deterministic system (this is a consequence of Proposition [l]) 
and reaches an e-neighborhood of h in finite time. Finally, the subpopulations which have a 
zero coordinate in h can be approximated by sub-critical branching processes. The durations of 
the first and third steps are proportional to log(A"), whereas that of the second step is bounded. 
The second inequality in (4.1) guarantees that, with high probability, the three steps of invasion 
are completed before a new mutation occurs. 


4.1.1. Invasion fitness. Given a population in a stable equilibrium that populates a certain set 
of traits, say M C X, the invasion fitness f(x, M ) is the growth rate of a population consisting 
of a single individual with trait x M in the presence of the equilibrium population h on M. 
In the case of the BPDL model, it is given by 

f(x, M) = b(x) - d(x) - ^ c(x, y)h y . (4.2) 

y&M 

Positive f(x, M ) implies that a mutant appearing with trait x from the equilibrium population 
on M has a positive probability (uniformly in I\) to grow to a population of size of order K ; 
negative invasion fitness implies that such a mutant population will die out with probability 
tending to one (as I\ —> oo) before this happens. 

Invasion fitness is a fundamental concept in the analysis of stochastic population models. 
We first generalize it to the case where fast phenotypic switches are present, for pure tumor 
evolutions, i.e. we ignore the T-cells and the cytokines in our model. 

Let us consider an initial population of genotype g (associated with l different phenotypes 
pi,... ,pe) which is able to mutate at rate to another genotype g', associated with k differ¬ 
ent phenotypes p\,...,p' k . The assumption ensures that no mutation occurs during the 
equilibration phase in the phenotypic space. 

Consider as initial condition n(0) = (n( 3jPl )(0),... >R( 9 , P£ )(0)) a stable fixed point, n, of the 
following system: 

( l l \ l 

h — di — ° ij n ( g , Pj ) — ^2 Sj d ) ^2 s A n (g,Pj)’ i = 1 ,... (4.3) 

3 =1 3 =1 / 3 = 1 

We write for simplicity bj = b(pi), d t = d(pi), Cij = c(pt,pj), Sij = s g (pi,pj), and later b[ = b(p , i ), 
d 'i = d( yP'i)i = c (p'vPj), Cij = cipiip'j), aC = Sgfip'^p'j). 

We want to analyze whether a single mutant appearing with a new genotype g' (and one of 
its possible phenotypes p(), has a positive probability to give rise to a population of size of order 
K. Using the same arguments as Champagnat et al. na i27], it is easy to show that as long 
as the mutant population has less than eK individuals (with e<l), the mutant population 

^By stable fixed point we mean that the eigenvalues of the Jacobian matrix of the system at the fixed point 
have all strictly negative real parts. 
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(r/ ,p\)..... (g',p' k ) is well approximated by a A;-type branching process, where each individual 
undergoes binary branching, death, or switch to another phenotype with the following rates: 


Pi -> PiPi 
Pi ->• P'j 


with rate 
with rate 
with rate 


b' 

d'i + £z=i 


for i,j G {1,.. ■, k}. 


(4.4) 


where (hi,..., n^) denotes the hxed point of (4.3). We 
the transition rates of an irreducible Markov chain on 
case where s- > 0, for all i, j G {1,..., k}. 

Multi-type branching processes have been analyzed 
Athreya and Ney [33]. Their behavior are classified in 


will assume that the switch rates ■ are 
{1,..., k}. The simplest example is the 

by Kesten and Stigum |3T, 32, [33] and 
terms of the matrix A, given by 


where 



( h 

^12 

4\ 




A = 

S 21 

h 


5 


(4.5) 


V4i 


fk / 




t 

dil 

1=1 

hi - 

k 

3=1 

i = 1, 

,k. 

(4.6) 


Note that fi would be the invasion fitness of phenotype i if there was no switch back from the 
other phenotypes to i (or if all switched cells would be killed). It is well-known that the multi¬ 
type process is super-critical, if and only if the largest eigenvalue, Ai = Ai(A), of the matrix A 
is strictly positive, meaning that if Ai > 0, the mutant population will grow (with rate Ai) to 
any desired population size before dying out. Thus Ai(A) is the appropriate generalization of 
the invasion fitness of a genotype: 

F(g',g):=\ 1 (A). (4.7) 

This notion can easily be generalized to the case when the initial condition is the equilibrium 
population of several coexisting genotypes. Note that this notion of invasion fitness of course 
reduces to the standard one of m if there is only one mutant phenotype. This settles the first 
step of our analysis, which is the invasion of the mutant. 


4.1.2. Towards a generalized Polymorphic Evolution Sequence. In fact, one can say more about 
how the mutant population grows. Write Zj l \t) for the number of individuals of phenotype pj 
existing at time t for this branching process when the first mutant is of phenotype pi. Then, 
for i , j < k, 

E (zf(t)) = [M(t)\ iJ (4.8) 

where M[t) is the k x ^-matrix 

M(t ) = exp(At). (4.9) 

Assume that the largest eigenvalue Ai(A) is simple and strictly positive. Let v and u be the 
left and right eigenvectors of A associated to Ai, normalized such that u ■ v = 1 and u ■ 1 = 1. 
The extinction probability vector q = (q\,... ,q k ) where qi = P (Z^\t) = 0 for some t) is the 
unique solution of the system of equations: 

l k ( i k \ 

d'i + ^2 CiWi + b[q. 2 + ^2 S ljs q 3 = qi d' + ^ cu n t +^ + 22 2 ) > i = l,...,k (4.10) 

i=i j=l V l=i 3=1 ) 

which has in general no analytical solution. Then the following limit theorem holds EH: 
hm (Z«(f),..., Z«(t)) e- A 'd = 


Wi • (vi,...,v k ) a.s. 


(4.11) 
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where (Wj)j=i ^ is random vector with non-negative entries such that 

P( Wi = 0) = Qi and E(Wj) = Uj. 


(4.12) 


In particular, conditionally on survival, the phenotypic distribution of the mutant populations 
converges almost surely to a deterministic quantity, which moreover does not depend on the 
phenotype of the initial mutant, namely 


V Z f® 

lim ——^—7TT- 

t—} OO r 7\ % ) 


e = i^ i; (t) e;=i- 


Vi,j' = 1,..., k. 


(4.13) 


For us, this implies the important fact that the phenotypic structure of the mutant population 


once it reaches the level eK > 0 is almost deterministic. Then, conditionally on survival, (4.11) 


implies that the time t £ k until the total mutant population reached eK is of order log(Ii). 
Moreover, the proportions of the k types of phenotypes converge to deterministic quantities 
given above, 


■ ■ ■, Z k (r £ K )) 


ev i 


ev k 


X,=ib £,=i v 3, 


, in distribution, as K —> oo. (4.14) 


Once the mutant population has reached the level eK , the behavior of the process can be 
approximated by the solution of the deterministic system: 

/ l l k \ l 

\a,Pi) = n (g,Pi ) I h — di - ^ s l3 — ^ c ij n (g,Pj) “ ^ ^i n (g',p') I + s A n (9>Pi)’ 

\ J=1 3 =1 J=1 ‘ / 1=1 

k k t 


n (9',p'i) n {g',p'i ) 1 ^ C \ 

1=1 1=1 


Cij n (g,Pi) I + s l* n (s'iPjO> (^-15) 
1=1 / 1=1 


'/h. </'•//,) 


(where i runs from 1 to t in the first line and from 1 to fc in the second one) with initial 
conditions in a small neighborhood of 


( n h?,Pi)(°)> ■ • ■. n O,pd(°)’ n (9Vi)(°)’ ■ • • > n (V,p()(°)) = Ui, • • •, 


£V\ 


£V k 


E k 

1=1 V 3 


E k 

1=1 ^i. 


( 4 - 16 ) 

If the system (4.15) is such that a neighborhood of (4.16) is attracted to the same stable fixed 


point, we are in the same situation as in Champagnat and Meleard m and get the generalization 
of the Polymorphic Evolution Sequence on the genotypic trait space. 


Observe that in the case when the system of equations (4.15) has multiple attractors, and 


different points near (4.16) lie in different basins of attraction, then for finite K, the choice 


of attractor the system approaches may be random. The characterization of the asymptotic 


behavior of the system (4.15) is needed to describe the final state of our stochastic process. 


This is in general a very difficult and complex problem, which is not doable analytically and 
requires numerical analysis. 

Figure [lO] shows examples where in a population consisting only of type (g,p) a mutation to 
genotype g' occurs. <j is associated to two possible phenotypes p\ and p 2 . 


Figure 10 (A) is realized with following parameters: 


bo = 5 
do = 0 
coo = 1 
coi = 1 

C02 = 0 


bi = 6 

d\ = 0 
cio = 1 
c n = 1 
C12 = 0 
n (9',pi)(°) =° 


b 2 = 6 
d 2 = 0 
C20 = 1 
C21 = 0 
C22 = 1 
n (9V 2 )(^ = ’ 


«12 = 0.1 
s 2 \ = 2 


(4.17) 


K = 200 
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A 



B 



FIGURE 10. Simulations for rare mutations in combination with fast switch¬ 
ing, where the number of individuals divided by 200 is plotted versus time. (A) 
The mutant phenotype P 2 has a negative initial growth rate but can switch to 
pi which has a positive one. The fitness of the genotype g' is positive. (B) The 
fitness of the mutant genotype g' is positive, although each phenotype has a 
negative initial growth rate. This is possible because an outgoing switch is a 
loss of a particle for a phenotype, but not for the whole genotype. 


In this case, p' 2 can switch to p\ but the back-switch is very weak, and we have f 2 < 0 and 
fi > 0 according to definition (4.6). The global fitness of the genotype g' is positive and close 


to /i. The coexistence of the two phenotypes depends on the cross-competition C 12 and C 21 . 


Figure 10 (B) shows an example of the case discussed above with the following parameters: 


6 0 = 5 
do = 0 
coo = 1 

coi = 1 
C02 = 0 

n (9iP) (°) = 5 


b\ = 6 
d\ = 0 
cio = 1 
cn = 1 
C 12 = 0 
n (9Vl)(0) = 1 


b 2 = 6 
d 2 = 0 
C 20 = 1 
C 21 = 0 
C 22 = 1 
n (9',^)(°) =° 


S 12 = 2 
S 21 = 2 


K = 200 


(4.18) 


For these parameters /1 and f 2 as defined in (4.6) are negative, but, due to the cooperation of the 


two phenotypes, the fitness of the genotype is positive and it invades with positive probability 


as indicated by the definition (4.7). Moreover, both phenotypes appear on a macroscopic level. 


4.2. Interplay of mutation and therapy. In the previous section we considered the proba¬ 
bility of invasion of a mutant when the resident population is at an equilibrium given by a stable 
fixed point. In the context of therapy, there are phases when populations shrink and regrow 
due to treatment and relapse phenomena. In the medical literature, there are frequent allusions 
to the possibility that such growth cycles may induce fixation of a “super-resistant mutant”, 
see e.g. [281 5] [29], It is important to understand whether and under what circumstances such 
effects may happen. 

Here we show an example where the appearance of a mutant genotype may indeed be en¬ 
hanced under treatment. We consider birth-reducing competition (BRC) between tumor cells. In 
such a case, a large population at equilibrium may encounter fewer births and hence mutations, 
than a smaller population growing towards its equilibrium size. 

Let us discuss in more detail how the birth-reducing competition can have a crucial effect on 
the mutation timescale. For the sake of simplicity we consider an example where the switching 
rates are set to 0. Consider a melanoma population ( g,p ) which is able to mutate to a fitter 
type of melanoma ( g',p')■ We allow for one T-cell population attacking the resident melanoma 
population since this is the simplest scenario where the effect of therapy in this context can be 
explained. As the presence of TNF-a only influences the switch between phenotypes, it does 
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not play any role in this example and we can set the corresponding parameters (£ w and d(w )) 
to zero. If /r A —> 0 as K —> oo the limiting deterministic system is given by: 


^(g,p) = n (.g,p ) i b (p) ~ d (p) ~ c b{p,p)n(g,p) ~ Cb(p,p')n { gf y) - t(z,p)n z ) 

\g', P ') = n (g',p') W) ~ d (p' ) “ c b(p',p')^(g\ P ') ~ Cb{p',p)n ig!p) ) 

h z = n z (b{z,p)n {g ^ p) - d{z )) (4.19) 

Note that the mutation term does not appear in the deterministic system and that the difference 
between birth-reducing competition and usual competition is not visible. The effects we are 
looking for are intrinsically stochastic and happen on time-scales that diverge with K. 

If the competition is only of birth-reducing type, then the total mutation rate of the popula¬ 
tion of type ( g,p ) at time t is 

™(^f (g,p)) ■= Pg 1 Kp) - c b {.p,p)^{g,p)\ + v? (g,p)K. (4.20) 


This is a positive and concave function of ( g,p ) on the interval [0, b(p)/cb(p,p)\, see Figure 
11 If the population is at equilibrium (without or before therapy), meaning (g,p) = h( 9)V ) = 
JJjyp) — d(p))/cb(p,p ), then the time until a mutation occurs is approximately exponentially 
distributed with parameter equal to 




( b(p) - Cb(p,p)n (gtP) ) n (3jP) = pfK ■ d{p) n (5iP) . 


(4.21) 


If d{p) is smaller than b(p)/2, then h( 9iP ) is bigger than b(p)/2cb(p,p) and m(h( 9p )) is not 
maximal. Smaller populations, more precisely in between d(p)/cb(p,p ) and h( 9lP ), have a higher 
total mutation rate. 

The interesting scaling of the mutation rate is K—y a > 0 as K —> oo. In this case, there 
is a number of mutations of order one while the population grows by 0{K) individuals. For 
lower mutation rates, no mutant can be expected before the population reaches its equilibrium, 
while for higher rates, mutations occur unrealistically fast. Since for /j. A -a 0 the mutation term 
does not appear in the deterministic system, the difference between BRC and usual competition 
is invisible. 


During therapy, a tumor which is close to equilibrium (similar to Figure 12 (B)) can shrink 


to a small size (similar to Figure 12 (A)): the introduction of T-cells in the system lowers the 


population size of melanoma, and the total mutation rate in the tumor population of type (g,p) 
can be larger during treatment, see Figure [l2| (C). This means that treatment could lead to 
earlier mutations and thereby accelerate the evolution towards more aggressive tumor variants. 



Figure 11. Initial total mutation rate of the population ( g,p ). 
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The simulations are obtained with the following parameters: 


b(p ) = 4 
dip) = 0.1 
Cb(p,p) = 3 
Cb(p,p') = 0.8 

n (9,P)(°) = L3 


b(p r ) = 6 
d{p') = 1 
Cb{p\p) = 0.8 
Cb(p ',£>') = 1 
n (ff',P')(°) = 0 


b(z,p ) = 20 
t(z,p) = 10 
d(z) = 6 

n z (0) = 0 or 0.1 


m((g,p ), = 1 

dg = 10- 3 

K = 10 3 


(4.22) 


Note that this example provides an interesting situation of interplay between therapy and 
mutation. By lowering the melanoma population, the T-cell therapy actually increases the 
probability for it to mutate to a potentially fitter and pathogenic genotype, which is not affected 
by the T-cells. 



(g,p) 


(g\p’) - 


- 

/ 

- 






Figure 12. Simulations of mutation events in a population, where com¬ 
petition is acting via birth reduction (parameters are given by ( 4.22| >). 
The number of individuals divided by 1000 is plotted versus time: Ef¬ 
fect for an initial population which is small (A), or at equilibrium (B) 
or under therapy (C) and (D). 


5. Discussion 

Therapy resistance is a major issue in the treatment of advanced stages of cancer. We have 
proposed a stochastic mathematical model that allows to simulate treatment scenarios and 
applied it to the specific case of immunotherapy of melanomas. Comparison to experimental 
data is so far promising. The models pose challenging new mathematical questions, in particular 
due to the interplay of fast phenotypic switches and rare driver mutations. First numerical 
results point to a significant effect of stochastic fluctuations in the success of therapies. More 
precise experimental data will be needed in the future to fit crucial model parameters. While 
our models describe the actions of individual cells and cytokines, they do not by far resolve the 
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full complexity of the biological system. In particular, they do not reflect the spatial structure 
of the tumour and its microenvironment. Also, the distinction of only two phenotypes of the 
tumour cells is a simplification. The same is true for the interaction with other immune cells 
and cytokines. This reflects on the one hand the limitation due to available experimental data, 
on the other hand the use of a model of reduced complexity also makes numerical computations 
and theoretical understanding of the key phenomena feasible. The rates entering as model 
parameters therefore have to be understood as effective parameters, e.g. the death rate of 
T-cells accounts for their natural death as well as the exhaustion phenomenon. In principle 
it is possible to increase the resolution of the model; this, however, increases the number of 
parameters that need to be determined experimentally which would pose a major challenge. 
Already at the present stage, the model parameters are not known well enough and are adjusted 
to reproduce the experimental findings. Some parameters that it would be very useful to see 
measured precisely are: 

• birth and death rates of tumour cells, both in differentiated and dedifferentiated types. 
Currently these are estimated from the growth rate of the tumour, but this yields only 
the difference of these rates; 

• killing rates of T-cells, both of the differentiated and the dedifferentiated tumour cells; 

• rates of phenotypic switches, both in the absence and the presence of TNF-ck; 

• death rates of T-cells and their expansion rates when interacting with tumour cells. 
Nevertheless, we see the proposed model as a promising tool to assist the development of 
improved treatment protocols. Simulations may guide the choice of experiments such that the 
number of necessary experiments can be reduced. The obvious strength of our approach is to 
model reciprocal interactions and phenotypic state transitions of tumour and immune cells in 
a heterogeneous and dynamic microenvironment in the context of therapeutic perturbations. 

The clinical importance of phenotypic coevolution in response to therapy has been recently 
documented in patients’ samples from melanomas acquiring resistance to MAPK inhibitors 
[35 j. Adaptive activation of bypass survival pathways in melanoma cells was accompanied by 
the induction of an exhausted phenotype of cytolytic T cells. This has important implications 
for the combinatorial use of cancer immunotherapy (checkpoint inhibitors like anti PD-1) with 
respect to scheduling. We envision that our mathematical approach will help to integrate such 
patient omics data with experimental findings to guide novel strategies. Of note and similar to 
our previous study, dedifferentiation of melanoma cells was identified as a major mechanism of 
escape from MAPK inhibitors (361 37], We recently dissected the molecular circuitries that con¬ 
trol melanoma cell states and showed how melanoma dedifferentiation governs the composition 
of the immune cell compartment through a cytokine-based crosstalk in the microenvironment 
|38[ 139]. Hence, malignant melanoma is a paradigm for a phenotypic heterogeneous tumour and 
a future goal is to incorporate this increasing knowledge of melanoma cell plasticity into our 
method to refine its capability to model complex interactions with immune cells. 

Importantly, phenotypic plasticity in response to therapy is a widespread phenomenon and 
non-small cell lung cancer (NSCLC) represents a prominent example. A subset of NSCLCs 
harbour activating mutations in the epidermal-growth factor receptor EGFR and respective 
small molecule inhibitors (EGFRi) are potent first line cancer drugs for this NSCLC subtype. 
However, tumours invariably relapse and genetic selection of subclones with the second site resis¬ 
tance mutation EGFRT790M is the major event. A substantial number of relapse tumours show 
remarkable transitions from an NSCLC adenocarcinoma to a neuroendocrine-related small cell 
lung cancer (SCLC) phenotype |40j . Given the recent success of immune checkpoint inhibitors 
in NSCLC }41j , it will be of clinical interest to investigate the phenotypic coevolution of immune 
cells in the context of NSCLC-SCLC lineage transitions. Again, our mathematical approach 
could represent a valuable tool to support this research. Finally, our results suggest that sto¬ 
chastic events play an unanticipated critical role in the dynamic evolution of tumours and the 
emergence of therapy resistance that requires further experimental and clinical investigation. 
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Appendix A. Pseudo-code 


The depth diagram of the algorithm we used to generate the simulations in this article is 
given in Figure 13 The pseudo-code is given in Algorithm ([Tj). 


We use the following notations: let D be some discrete set and X a "D-valued random variable, 
then X sampled according to the weights {w(i),i G D} means that P(X = i) = w(i)/ Yliev w i})- 



Figure 13. Depth diagram of the program. 
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Data: initial conditions: Vq £ M K (T), number of iterations: Alterations, parameters described in Section 1 
T 0 < 0, t— i>o , k < 0 

While k < Alterations do 


for x £ Supp{yT k ) do 

if x = (g,p) egxV then 

K 


B(x) £- Ku Tk (g,p) 


KP) ~ E(g,p )e supp(^ ) c b(P,P)v Tk (g,P) 


D(z) <- Kv$ k {g,p) (d(p) + b(p) - E( 9 ,p )e supp(^) c b {p,p)v£ k (g,p) 


+ T,^p)eSu PP (^ k )<P’P>T k {g,p)j, 

T(x) <r- Kv^ k {g,p) E 

zGSupp(i/^) t{z,p)v% k {z), 

S(x) <— Kux k (g,p) T,pev (s 3 {p,P) + E„, € su PP (**) 

P(x) «— 0, 

if x = z £ 2 then 

fl(a;) <-jr«# fc (*)6(z), D(x)^Kv* k (z)d(z), T(x) <— 0, 5(®) «- 0 

P(x) £- AVrJz) E (s , P )6 Su pp(^|f) b( z >p) l 'T k {giP), 


if .x = w £ W then 

[ B{jc) t— 0, _D(x) <—AV|f (w) d(w), T(x) <— 0, S^x) t— 0, P(x) t—0, 
TotalTraitRate(x) <— B(x) + D(x) + T(x) + P(x) + S(x) 

TotalEventRate t— E^gSuppO^ )TotalTraitRate(x) 

^k ' 

Sample t ~ £xp(TotalEventRate) 

Tk+1 t - + t 

Sample x c hosen £ T according to the weights {TotalTraitRate(x), x £ Supp(i^f ))}. 

Sample the event E £ {Birth, Death, Therapy, Production, Switch} according to the weights 

{.£? (Xchosen ) , -D(Xchosen), T(Xchosen), -P(Xchosen), S (.Xchosen ) } 

case E = Birth 

if X chosen — (g,p) £ Q x V then 

Sample B £ {No Mutation, Mutation} according to {1 — p g , p g } 

case B = No Mutation 

L U T k + 1 ^ V T k + iM^chosen 

case B = Mutation 

Sample (g,p) according to m((g,p), ( g,p )) 

V T k + 1 <r- VT k + ^<5(9,p) 


else 

L U T k+ l < ^ U T k + chose, 


case E = Death 

i t— ZX? — b&x V, 

[_ 1 k+ 1 J-k K ^Chosen 

case E = Therapy 

Note that x c hosen = (g,p) for some (g,p) £ Q x V in this case. 

Sample z £ Z according to the weights { t{z,p)vx k (•?), z £ Supp(zx}f ) n Z } 

u r k+1 t- »T k - Tc&( g ,v) + E u , 6 w^™ U ( 2 : ’P)^' 5to 

case E = Production 

Sample (g,p) £ Supp(^) according to the weights {&(x c hosen,p)z'|^ {g,p), (g,p) £ Supp(z/|l)} 

u T k+ i ^T fc + K ^chosen + EwgW (Xchosen, p) 

case E = Switch 

Note that x c hosen = (<?,p) for some (<?,p) £ Q x V in this case. 

Sample p £ V according to | s 3 (p,p) + E^gsupp^) (P,P)^T k (w),p £ p| 

^T fc + 1 ^ u T k - ~k S (g,p) + -R 5 (g,p) 

k = k + 1 


Algorithm 1: Pseudo-code of the Gillespie algorithm used for generating the figures. 
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